function [x, y] = GPStoXY(lat, lon, ref_lat, ref_lon)
	% input GPS and Reference GPS in degrees
    % output XY in meters (m) X:North Y:East
    CONSTANTS_RADIUS_OF_EARTH = 6371000.;
    
    lat_rad = lat*pi/180;
    lon_rad = lon*pi/180;
    ref_lat_rad = ref_lat*pi/180;
    ref_lon_rad = ref_lon*pi/180;

    sin_lat = sin(lat_rad);
    cos_lat = cos(lat_rad);
    ref_sin_lat = sin(ref_lat_rad);
    ref_cos_lat = cos(ref_lat_rad);

    cos_d_lon = cos(lon_rad - ref_lon_rad);

    arg = min(max(ref_sin_lat * sin_lat + ref_cos_lat * cos_lat * cos_d_lon,-1.0), 1.0);
    c = acos(arg);

    k = 1.0;
    if abs(c) > 0
        k = (c / sin(c));
    end

    x = k * (ref_cos_lat * sin_lat - ref_sin_lat * cos_lat * cos_d_lon) * CONSTANTS_RADIUS_OF_EARTH;
    y = k * cos_lat * sin(lon_rad - ref_lon_rad) * CONSTANTS_RADIUS_OF_EARTH;